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ABSTRACT 


We present calculations of AGN winds at ~ parsec scales, along with the associated 
obscuration. We take into account the pressure of infrared radiation on dust grains and the 
interaction of X-rays from a central black hole with hot and cold plasma. Infrared radiation 
(IR) is incorporated in radiation-hydrodynamic simulations adopting the flux-limited diffu- 
sion approximation. We find that in the range of X-ray luminosities L=0.05 - 0.6L e dd, the 
Compton-thick part of the flow (aka torus) has an opening angle of approximately 72° — 75° 
regardless of the luminosity. At L > 0.1 the outflowing dusty wind provides the obscuration 
with IR pressure playing a major role. The global flow consists of two phases: the cold flow 
at inclinations 9 > 70° and a hot, ionized wind of lower density at lower inclinations. The 
dynamical pressure of the hot wind is important in shaping the denser IR supported flow. At 
luminosities <0.1Z/ e( jd episodes of outflow are followed by extended periods when the wind 
switches to slow accretion. 


1. Introduction 

Active galactic nuclei (AGNs) are among the most fascinating objects in the universe: 
various energetic processes interplay in a volume a few cubic parsecs around a supermassive 
black hole (BH); radiation is likely decisive in connecting the accretion and feedback modes 
of AGN. Dusty, radiation driven flows are very effective in removing large masses of accreting 
gas, returning it back the galaxy, and naturally limiting the growth rate of a supermassive BH 
by the feedback with the host galaxy. Radiation also provides a non-linear causal connection 
between distant parts of an AGN itself, coupling together vastly different physical scales. 
Despite the pivotal role of radiation in the first principle modeling of AGN and the significant 
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progress in theoretical modeling over the last 10 years, proper incorporation of radiation into 
numerical models still presents a significant computational challenge. 

Well established phenomenology places AGNs into two distinct classes which roughly 
differ by the presence (type I) or absence (type II) of the broad emission lines in the optical 
and UV spectra. According to AGN unification paradigm (Rowan- Robinson 1977; Antonucci 
& Miller 1985; Urry & Padovani 1995), this dichotomy is an artifact of a pure geometrical 
origin. In order to be observed as type II object, an AGN must be viewed at high inclination, 
in extreme cases close to the equatorial plane, while in the case of a type I AGN the observer’s 
line of sight is at lower inclination, i.e. closer to the symmetry axis. The presence of 
geometrically thick toroidal obscuration is pivotal to the unification scheme of AGNs. 

The discovery of broad permitted emission lines in the polarized spectrum of the nearby 
Seyfert II galaxy NGC 1068 Antonucci (1984), and Antonucci & Miller (1985), provides 
evidence that a type I nucleus is harbored inside the obscuring ring of matter. The observed 
polarization is likely due to scattering of light emitted by the type-I core and then being 
scattered within an outflow emanating from the same nucleus. 

Observations of AGNs and quasars reveal that dust in large quantities can be present 
within a few parsecs from the BH, and a significant part of X-rays and UV are reprocessed 
into infrared (IR). For example, in luminous quasars almost half of the bolometric luminos- 
ity can be reprocessed into IR. This results in a 1-100/im ’’hump” in the quasar spectral 
energy distribution (SED) Richards et al. (2006). In nearby AGNs multi-temperature dust 
distributions are observed directly by the Very Large Telescope Interferometer (VLTI). Such 
mid-infrared observations of the prototypical Seyfert II galaxy NGC 1068 (Jaffe et al. 2004) 
and the Circinus galaxy (Tristram et al. 2007) depicts a picture of a dusty, clumpy plasma 
exposed to external illumination. 

Various models have been proposed to explain the real or apparent thickness 
of the torus: 

1) a ” warped disk model”, these are global warps in a locally geometrically 
thin accretion disk (i.e., Phinney 1989; Sanders et al. 1989), 

2) i) models of a quasi-static, rotating torus involving ’’clouds” (i.e., Krolik & 
Begelman 1988; Beckert & Duschl 2004; Nenkova et al. 2008). Vertical thickness 
of the torus is supported via large velocity dispersion of clouds, and small-scale 
magnetic field is required to maintain elasticity of clouds during collisions; ii) 
models which rely on infrared radiation pressure on dust to keep rotating torus 
vertically thick (Krolik 2007). Most of IR radiation in this models comes from 
reprocessing of external X-ray, and UV radiation. 
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3) a ’’wind torus” A collar of outflowing matter may provide obscuration, 
i) MHD wind: possibility is that a global magnetic field can be involved in the 
form of an MHD wind (Konigl & Kartje 1994; Elitzur & Shlosman 2006) or di- 
rectly supporting a quasi-static torus (Lovelace et al. 1998), or in a combination 
of radiative and magnetic driving (Emmering et al. 1992; Everett 2005; Keat- 
ing et al. 2012); note that sufficiently strong (of the order of the equipartition 
one) poloidal magnetic field can also help to solve the problem of how the an- 
gular momentum is removed from accreting matter, however the origin of such 
highly ordered field remains unclear ii) an outflow (or failed wind) driven by IR 
radiation pressure on dust (Dorodnitsyn et al. 2011, 2012). All these models 
have their benefits and drawbacks. Models investigating the properties of a 
self-gravitating, non-radiative, dissipationless torus through N-body simulations 
have been made by Bannikova et al. (2012). 

The accreting dusty gas at ^parsec distances from the BH is susceptible to 
self-gravitating instability; since it is exposed to intense radiation strong gravito- 
thermal instabilities can develop. 

In the outflowing wind, the non-linear phase of thermal instability is associ- 
ated with ’’clumps” (i.e., Elitzur 2008). Presumably such clouds are confined by 
small scale magnetic field. Whether they are secondary to the flow and behave 
more like ’’waves”, i.e. temporary structures evolving on a dynamical or radia- 
tion times scales, or real long-lasting ’’clumps” remains unresolved. All models 
involving clouds need to explain how to avoid too much dissipation through the 
cloud-cloud collisions, which would otherwise raise the temperature to a fraction 
of the virial temperature T v i r g . In the present studies we do not consider the 
effect of magnetic field (both ordered and chaotic). 

It has long been recognized that the geometrical thickness of the obscuration 
presents a significant problem. In order to explain why the torus is geometrically 
thick, models have been constructed which rely on gas pressure, magnetic forces, 
and on radiation pressure either alone or in combination with other mechanisms. 
If one relies on gas pressure to support vertical thickness of the torus then virial 
arguments suggest that the temperature of such gas should be of the order of 
T v [ r g = 2.6 x 10 6 M7/r pc K, where M 7 is the black hole (BH) mass in 10 7 M o , and r pc 
is the distance in parsecs. Such high temperatures are not compatible with the 
existence of dust. 

As mentioned before, the presence of a small scale magnetic field can provide 
necessary elasticity to cloud-cloud collisions. However, this argument is not 
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without a flaw. It is natural to assume that turbulence develops. If so, this leads 
to the magnetic field being of the order of the equipartition field. A result is 
that the dissipation of magneto-sonic waves will raise the electron temperature 
to be of the order of T vir g . 

The connection between a toroidal obscuration, UV absorption and warm 
absorbers can be self-consistently assessed by global multidimensional simula- 
tions. The prequisite for such an endeavor is to include various physical effects 
and couple them to radiation. Additionally why such calculations are currently 
beyond the reach is that high numerical resolution is needed to resolve multi- 
phase behavior of the outflowing gas. These are the reasons why different works 
including present studies assess smaller aspects of the bigger picture. 

The connection between a Broad Line Region and dusty outflow was con- 
sidered by Czerny & Hryniewicz (2011). They suggested that the gas of the 
dusty accretion disk is weekly coupled to the BH due to radiation pressure on 
dust resulting in a dusty wind which gives the low- ionization part of the broad 
line region. Bottorff et al. (2000) considered the relation between the observed 
properties of UV and X-ray absorbers with those predicted from a magneto- 
hydrodynamical self-similar wind model. They conclude that in case of the 
well-studied Seyfert 1 galaxy NGC 5548 the intrinsically clumped absorber is 
biased toward smaller radii in case of warm absorber and are likely located at 
higher distances in case of UV absorption features. 

Other models are based on the possibility that the pressure of the infrared radiation on 
dust may itself be enough to support the vertical thickness of the torus. The idea of a static 
infrared support was explored in Krolik (2007) where an approximate analytical model for 
the torus was constructed showing that there is enough radiation force to support a static 
vertically thick torus. However, it was argued in Dorodnitsyn et al. (2011) (Paper I), that it is 
very unlikely that the IR supported torus would be static. It was shown that the equilibrium 
between rotational, gravitational and radiation forces cannot be maintained if the temper- 
ature in the toms exceeds that of T v ir,r — 312 — 987 (n/10 7 Mjr^) 1 ^ 4 K, 

where n is the number density, necessarily resulting in the torus being an outflow. In the 
following paper (Dorodnitsyn et al. 2012, Paper II) we supported this idea with radiation- 
hydrodynamics simulations of the outflowing accretion disk wind, driven by the radiation 
pressure on dust. The distribution of IR radiation was calculated using a flux-limited diffu- 
sion approximation. 

An important simplification was made in Papers I and II: Only the infrared-driven 
part of the wind was considered. The argument was that this portion of the wind is far 
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enough from the BH so all the conversion of X-rays to IR already happened outside the 
computational domain. That allowed us to consider only one group of radiation (the IR). 
The radiation temperature at the boundary was calculated to match the input energetics 
of soft X-rays. In other words we did not treat explicitly the reprocessing of soft X-rays 
into IR. Correspondingly, we were not able to follow the interaction of the X-ray heated hot 
component of the flow with the cold, dusty IR supported component. 

In this paper we calculate the wind structure and explicitly treat the interaction between 
X-rays and gas with subsequent conversion of X-rays to IR. This allows us to follow in 
dynamics the transition from hot, photoionized to cold, dusty and IR supported flow. To 
make calculations numerically tractable we made several major simplifications in treatment 
of the transfer of X-rays as well as their interaction with the gas. For example, transfer 
of X-rays is reduced to simple attenuation with no multiple scatterings taken into account, 
and we approximately calculate the heating of dust grains by X-rays. The rest of our 
computational radiation-hydrodynamics framework remains unchanged from Paper II: the 
transfer of the infrared radiation and its interaction with matter is treated in a flux-limited 
diffusion approximation. The calculations are time-dependent and 2.5D, i.e. 3 dimensions 
in cylindrical coordinates and assuming axial symmetry. 

The plan of this paper is as follows: After this introduction the flow of this paper goes 
through Section 2 where we review basic assumptions of our model, including the description 
of the radiating fluid of the AGN dusty wind by means of radiation hydrodynamics equations. 
Interaction of the cold and hot components of the wind with X-rays are described in Section 
3 including corresponding approximations. In Section 4 we describe the numerical setup, 
initial and boundary conditions. The results are summarized in Section 5 while the final 
discussion is postponed to Conclusions, where we highlight prospects and limitations of our 
approach to the problem of AGN winds. 


2. Assumptions and physical model 

The flow in which the continuum radiation plays an important role is described by the 
system of radiation hydrodynamics equations. To first order in v/c, these can be formulated 
in the following form (Mihalas & Mihalas 1984): 
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where quantities related to matter: p , p, and e are the material mass density, gas pressure 
and gas energy density, v is the velocity; quantities related to radiation are the frequency- 
integrated moments (for definitions see Appendix A): E is the radiation energy density, F 
is the radiation flux, and P is the radiation pressure tensor; xpj xe are the Planck mean 
and energy mean absorption opacities (in cm -1 ), c is the speed of light, B = aT 4 / n is the 
Planck function, and a = ac/A is the Stefan-Boltzmann constant, T is the gas temperature; 
other notation include the convective derivative, D t = A + v • V, and Vv : P denotes the 
contraction ( djVi)P l A Notice that all the dependent variables in (l)-(4) are evaluated in the 
co-moving frame. Notice that E = aT j 4 , where T r is the radiation temperature. Interaction 
between X-rays and the gas is incorporated in the last term of the gas energy equation (3), 
where T and A are respectively heating and cooling rates. The equation of state for the gas 
is assumed to be polytropic: p = Kpx, where 7 = 1 + 1 /n, and n is the polytrope index, 
and p = (7 — l)e. The gas is one-component, one-temperature T = p/j/ plZ, where p is the 
mean molecular weight per particle, 7 Z = 8.31 x 10 7 ergK -1 g -1 is the universal gas constant 
and plasma with 7 = 5/3 is assumed to constitute the flow. Three components of the flow 
velocity vr, v z , and v ( t, are calculated in cylindrical coordinates z, R, assuming azimuthal 
symmetry, d ( p — 0 . 

Frequency-independent moments E, F, which appear in the above set of RHD equations 
are obtained by calculating angular moments from the frequency-integrated specific intensity, 
I(r,Q,u,t) and are given in Appendix A. Forces which are explicitly taken into account 
include radiation, g raf j pressure and gas pressure, p -1 Vp. 


These equations should be supplemented by the radiation force, g ra( j, which is calculated 
from the following relation 


Srad 


1 Xf F 
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where Xf = Xa + Xt is the total flux mean opacity consisting of absorption opacity, Xa an d 
the Thomson scattering opacity, xt- As in Papers I, and II we did not differentiate between 


7 


Xf, Xp and X'e- Radiation pressure on dust is the most important from the dynamical point 
of view (Paper I). Notice that in the infrared, the Rosseland mean opacity, k c ] = Xd/ P of 
dust with the temperature 10 2 — 10 3 K, is approximately 10 — 30 times larger than that 
of the electron Thomson opacities (Semenov et al. 2003). In this paper when calculating 
radiation pressure force in (5) we take into account only IR pressure on dust. Including 
UV force (in lines and on dust) is left for future studies. We adopt a simple procedure for 
the calculation of g ra< p if the temperature of the dust, T <: \ > T su b, where T su b is the dust 
sublimation temperature (see further in the text) the force is assumed to be zero. 

Our description takes into account two types of radiation: infrared and X-ray. 

The transfer of X-rays is treated in a single stream approximation. This implies that 
X-rays are simply traced from the point source (i.e. corona) adopting e~ Tx attenuation, 
where 7x is an optical depth with respect to X-ray absorption (Dorodnitsyn et al. 2008). 
The distribution of the infrared radiation is calculated in a Flux Limited Diffusion approxi- 
mation. Detailed discussion about the validity of equations (l)-(4) is given in Paper II. We 
briefly remind that the closure relation between F and E is obtained adopting the diffusion 
approximation: 


F = -DVE, (6) 

where if optical depth is large r> 1, the diffusion coefficient reads: 

D = cA, (7) 

where A = 1 /{up) is the photon mean free path, and k = x! P- From (5), (6), and (7) one 
can see that in the diffusion regime the radiation force does not explicitly depend on the 
opacity. 

The diffusion approximation tacitly assumes that optical depth r > 1. Most of the 
torus where IR pressure is important indeed has > 1 where r c ] is the optical depth of the 
gas-dust mixture in the infrared. 

If r < 1 the diffusion approximation should be modified so a correct limiting behavior 
at r < 1 is regained. One can see that without such modification when r<l the mean free 
path, A — » oo, and D — » oo, and |F| — » oo which is in contradiction with a free- streaming 
limit where it should be |F| — >■ cE (7). That is, when optical depth becomes small, or when 
p — » 0, the standard diffusion approximation is no longer applicable. 

In order to describe correctly regions of small r, the standard approach is to adopt 
the flux-limited diffusion approximation (Alme & Wilson 1974; Minerbo 1978; Levermore 


& Pomraning 1981). In this approximation A is replaced by A* = A A, where A is the 
flux limiter. The flux limiter we adopted in Paper I, II and in the current work is that of 
Levermore & Pomraning (1981): 


= 2 + ^LP /on 

6 + 3_R lp + i?L P ' J 

where ATp = A \WE\/ E. One can see, that if r — » 0, then i? L p — > oo, and |F| ~ cE; and if 
r 1 i?LP — > 0 and A — » 1/3. 


3. Interaction of radiation and the wind 
3.1. Photoionization equilibrium 

The premise of the current work is that the deposition of energy at parsec scale is 
dominated by external sources, i.e. by UV and X-rays from the central black hole and inner 
accretion flow. Due to the high opacity of the dust-gas mixture, a significant part of the 
UV and soft X-rays is absorbed and reprocessed into infrared in a thin layer of thickness, 
81 /Ripe — 1.3 x KT 3 nf 1 (Paper I). 

The dynamical time within the flow is usually much larger than the characteristic time 
of the photoionization and recombination. Thus, the ionization balance is determined by the 
condition of photo-ionization equilibrium. The condition of a photo- ionization equilibrium 
in the wind implies that the state of the gas can be be reasonably parameterized in terms 
of one ’’ionization parameter”, £ -the ratio of radiation energy density to baryon density 
(Tarter et al. 1969). In a popular form it reads 


£ = 4 7r F x / n ~ 4 • 10 2 • / x T M 6 / (N 2 3 r pc ), (9) 

where F x is the local X-ray flux, L x is the X-ray luminosity of the nucleus, and n is the 
gas number density; V 2 3 is the column density in 10 23 cm -2 , / x is a fraction of the total 
accretion luminosity Lbh available in X-rays, T is a fraction of the total Eddington luminosity 
L e dd = 1.25-10 44 Mq, where M 6 is the mass of a BH in 10 6 M o . We adopt a simple attenuation 
model in which the local deposition of energy is proportional to F x = L x e _Tx /(47rr 2 ), where 
the optical depth at soft X-rays is integrated over a straight line connecting the point source 

corona with a given parcel of gas r x = / K x p dl , and k, x is the X-ray opacity. 
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3.2. Hot and cold gas 

To calculate photoionization balance, we take into account Compton and photo-ionization 
heating and Compton, radiative recombination, bremsstrahlung and line cooling. The rates 
of these processes are calculated making use of the XSTAR photo-ionization code (Kallman 
& Bautista 2001) for the incident spectrum which is a power law with energy index a. As the 
direct incorporation of the XSTAR subroutines into the hydrodynamics code is not feasible, 
the key ingredient of our method is to make use of approximate analytic formulas for the 
heating and cooling rates. These formulae approximate extensive tables of heating-cooling 
rates obtained from XSTAR; these approximations are similar to those of Blondin (1994), 
although modified by Dorodnitsyn et al. (2008) to incorporate better atomic data and power 
law ionizing continuum with a— 1: 


Tic (erg cm 3 s -1 ) = 8.9 • 10 -36 £ (T x — 4T), (10) 

for the Compton heating - cooling; 

T x (erg cm 3 s” 1 ) = 1.5 • 1(T 21 f 1/4 T~ 1/2 (T X - T)T~\ (11) 

for the photo-ionization heating-recombination cooling, and for the bremsstrahlung and line 
cooling: 


A(erg cm 3 s -1 ) = 3.3 • 10“ 27 T 1/2 

+ (4.6 • 10“ 17 exp(— 1.3 • l0 5 /T)^ ( “ a8 “°- 98 “ ) T-l/2 + 10“ 24 ) 5. (12) 

Formulae (10)-(12) were found to be in reasonable (25%) agreement with numerical results. 

A variety of physical and chemical processes influence the transformation of X-rays 
into IR if dense dusty, molecular gas is illuminated by external X-rays (Maloney et al. 
1996; Krolik & Lepp 1989). Rather than attempt to incorporate all these processes in our 
radiation hydrodynamics framework, in the current paper we adopt several approximations 
which reproduce the qualitative behavior of low temperature gas exposed to external X- 
ray radiation. At column densities, N smaller than 10 24 cm -2 photons with energies below 
10 keV interact with the gas through photoelectric absorption while Compton scattering 
dominates at higher energies. The total ionization per particle at a given point in the 
wind is proportional to £ e ff, the ionization parameter which takes into account the frequency 
dependent attenuation of X-ray flux. It is shown by Maloney et al. (1996) that heating 
and cooling rates of such X-ray illuminated molecular gas depend on the modified ionization 
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parameter: £ e g = = 1-26 x 10 _1 F x /(n 8 A^ 22 9 ), and several regimes are approximately 

distinguished: 1) The regime of highly ionized gas £ e ff ^>1. In this regime for the heating 
and cooling rates we make use of approximations (10) - (12); 2) If 10 -3 -C £ e ff 1 the gas 
is primarily atomic, and we also adopt heating-cooling rates (10) - (12); 3) If the ionization 
parameter, is small (£ e fr -C £ m = 10 -3 ) the gas is largely molecular. This is where the most 
complexities due to chemistry and radiative interaction with the dust occur. The proper 
incorporation of this regime calls for consideration of all the radiative, chemical and dust 
network of processes. This is challenging even without coupling to hydrodynamics. 

In the regime when £ e ff < Cm, we assume that the molecular gas is in radiative equi- 
librium with dust, i.e. T g = T d where equilibrium dust temperature is found from the 
attenuated X-ray flux, T x . Thus we approximately take into account that dust can directly 
reprocess X-rays to IR contributing to the IR energy density, E. We approximately take this 
contribution into account, using equation (see, Appendix B): 

^ =n d ca d (T d - T r ) , (13) 

Dust temperature is found from the approximate relation, T d = (4F X / (ac)) 1 / 4 (notice that 
E = aT 4 ), where Fx is the local attenuated X-ray flux. If the dust temperature T d > T su b, 
where T su b is the dust sublimation temperature then it is assumed that dust is destroyed 
and no conversion to IR occurs. 

In our approach, the main difficulty in calculating the energy density of IR radiation 
from X-rays resulting from a single fluid approximation. By having only one species (gas) 
and not evolving a separate dust component, we can mimic the interaction between dust 
and X-rays, and calculate Td by means of equation (13), while additionally assuming the 
relation between Td, and T in a dusty- molecular phase, in which, for simplicity we assume 
T = T r , with the later calculated from(13). Other dust characteristics in (13) are as follows: 
n d = pfd/ m d is the dust grain number density, where f d is the dust-to-gas mass ratio, and m d 
is the dust grain mass, and a d = Tir\ is the dust cross-section, and r d is the dust grain radius. 
Throughout our calculations, we adopt the following fiducial values for the parameters of 
dust: typical dust grain radius, r d = 1 x 10 -4 cm, the density of the dust grain pd = 1 g cm -3 , 
fd = 1 x 10 -2 , and above the sublimation temperature T su b = 1500 K the dust is assumed to 
be destroyed. Equation (13) for T d (Fx ) approximately describes how X-rays are converted 
into IR through reprocessing by dust. 

If seff T bin and the T < T su b we assume that the exchange between T r ( E ) and T(e) is 
described by the balance between D t (E/p) and the last two terms in equation (4). The gas 
temperature T can contribute to IR and vice versa. For example, gas can contribute to the 
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X-ray-IR conversion through cooling. Cooling can be due to radiative processes, i.e. the low 
energy limit of (11)-(12) and due to ’’adiabatic” cooling due to fluid motion. 

The update of E from (13) is performed adopting operator splitting in the source step, 
i.e., before the update of the equation (4) is done. During the update over dt, the new 
E r (t + dt) is calculated from (13). 


4. Methods 

To solve the system of equations (l)-(4) we expanded the radiation hydrodynamic frame- 
work described in our previous papers. The additions include the interaction of X-rays with 
dust grains and with gas. The coupling between the X-ray heating and cooling with hydro- 
dynamics has been implemented and tested in our earlier works (Dorodnitsyn et al. 2008). 
Here we have added physics which includes the direct interaction of X-rays with dust grains. 

The radiation- hydrodynamic equations are solved using the ZEUS code (Dorodnit- 
syn et al. 2012) with extensions taken from Dorodnitsyn et al. (2008). The radiation- 
hydrodynamic part is designed and tested in (Dorodnitsyn et al. 2011, 2012) and conforms 
with methods and structure of the original ZEUS code (Stone & Norman 1992). 

In the following we briefly review the organizational structure of the code. Similarly 
to the way pure hydrodynamics is treated in ZEUS, the radiation part is also split into the 
source and transport steps. Artificial viscosity is adopted to resolve discontinuities in the 
flow. In the source step the difference equations for dE/dt and de/dt are solved, and in the 
transport step the advection of E and e is calculated. The update of E and e is calculated 
adopting fully implicit scheme (see Paper II); this includes the solution of the diffusion 
problem and the e E update i.e. the update of gas and radiation energy densities due to 
radiation-gas interaction. 

Similarly, in the X-ray-matter interaction step, heating y x — » e, and cooling e — > 7 X 
due to Compton, radiative recombination, bremsstrahlung and line cooling is also done fully 
implicitly in a separate update. 

We make use of the R, z cylindrical grid, which extends from R m to R out in radial and 
from Z\ n = 0 to z out = R 0 ut in the vertical direction. The lower boundary at z- m = 0 is 
intended to represent a thin accretion disk, which serves as an infinite reservoir of gas and 
dust. Wind can escape from all other boundaries, in principle, though the centrifugal barrier 
effectively prevents escape from the inner boundary. A fixed grid of size 200 x 200 is adopted 
for all calculations. Hydrodynamic boundary conditions (BC) are as follows: at the left 
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(innermost), right, and upper boundaries we adopt outflowing boundary conditions. 

At the equator we fix the inflow speed v 0 j at a small fraction of the speed of sound leaving 
the density to adjust according to the continuity equation. Hereafter, index j varies over R 
and index i over z. We allow pi S j to vary and fix v is j , that is p ls j = Vi S+ ijPi S+ ij/vi S , and = 
const , provided v is <C v ls+ 1;S , where v ls+ \ jS is the sound speed; everywhere in the simulations 
we assumed v is /v is+ i iS = 0.01. Notice that this is different from the implementation of the 
equatorial BC in Paper II where it was assumed p ls ,j ~ R~ p and velocity was found from the 
continuity equation and only checked after the fact that <C Uj S)S . Present implementation 
of boundary conditions allows the disk to be flexible choosing mass loading of the wind. The 
initial distribution of density at the equator was taken to be Por^ 1 , where p 0 = fn p n 0 is the 
characteristic density, m p is a proton mass. We find that numerical solution quickly forgets 
about initial equatorial distribution of p, and that solutions are not sensitive to the choice of 
p. Theoretically, since pi S j is allowed to change, the final solution should only weakly depend 
on initial equatorial distribution of p. Practically, however, this is true only for those parts 
of the disk where the wind is developed. We also find that the present implementation of BC 
produces a more stable solution near the equator, thus somewhat relaxing the restriction on 
the hydrodynamical time step. 

The radiative BC are free-streaming F ~ cE at all boundaries except at the equator 
where zero flux BC are implemented. 


5. Results 

5.1. Dynamics of the wind 

Deposition of radiation energy into the flow is controlled by the parameter t] = L\ M \ / Ly <: \d . 
We calculate models for p = {0.05,0.1,0.3,0.4,0.5,0.6}. The innermost boundary of 
the computational domain is located at radius, Ro = 0.5 pc. Other parameters are 

A/rh = 1 x 10 5 * 7 Mq, /x = 0.5, and the initial scale density n 0 = 1 x 10 8 * * cm -3 . The evolution 
of the radiative flow was followed for 10-30 dynamical times, to ~ 1.5 x 10 11 rp/ 2 Mf 1/2 s. 

In Papers I and II the temperature of the IR radiation was prescribed at the boundary, in 

a physically motivated way (to match the energetics of X-rays, / x Tboi)- The results obtained 
in Paper II are generally in a good agreement with the present modeling. However, explicitly 
including X-rays into calculations reveals several important features of the flow which were 
absent in previous studies. 

One of the results of the present study is that the pressure of the hot photo-ionized 
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component is very important in shaping the IR flow. Excitation of such a hot, photoionized 
wind automatically mirrors the creation of the IR flow. The interface between the two 
components of the flow is seen in all simulations presented in this paper. In all cases it 
takes the form of oblique shocks or contact discontinuities. A supersonically moving hot gas 
obliquely hits the slowly moving IR flow creating a discontinuity (in most cases an oblique 
shock wave). 


5.1.1. Model with L = 0.6L e dd 


Figure 1 shows the development of the radiation-driven disk wind for T = 0.6. The 
outflow develops quickly after about one dynamical time, f 0 — 4 x 10 4 yr. The mass-loss 
rate is (M) ~ 0.1 — O.2M 0 yr _1 on average, and peaks at (M) ~ lM 0 yr _1 in the period of 
time between 4 x 10 3 yr and 5 x 10 3 yr. The total computation spans from 0 to 8 x 10 3 yr. 
These averaged peak numbers include both failed wind and wind which has enough energy 
to escape to infinity. 


The energetics of the wind is measured by computing the kinetic output of the wind 
at the outer boundary of the computational domain E: i>km — 2Li tni /A/, where L^ in = 

J pv 3 /2dT. is the kinetic luminosity of the wind. For a radiation-driven wind we expect low 

values of Akin/Aboi! h ere we have Lkm/Tboi — 5 x 10 -5 — 2 x 10 -4 . 

The averaged dynamics of the wind, is described by the the average bulk velocity of the 

flow (v) = / pvdV / / pdV, where V is the total or partial volume occupied by the flow. 
Jv 


iv 


Here, the dense part of the flow has (v z ) — a few x 1km s 1 , and (vr) ~ a few x lOkms 1 . 
The maximum velocity reaches lOOOkms -1 in R, and over 600kms _1 in the z direction. 


Figure 2 shows a color plot of the distribution of the gas temperature. The gas spans 
a vast range of temperatures from being in a cold molecular state to photo- ionized, high 
temperature warm absorber flow. An apparent feature of our models is that the hot wind 
consists of large scale inhomogeneities. The cold flow doesn’t show such large scale structure. 
Notice that we have a very simple test if the gas is in the cold, molecular-dusty phase, 
£eff < £m, and that if hot component is not shown the density structure is much more 
pronounced in the temperature plot. 


It is of interest to address the question of whether there can be the cold 
gas in a hot wind as well. However, the limited resolution of our studies does 
not allow us to provide a reliable answer to this question. We do not find the 
coexistence of hot and cold phases of gas in any appreciable quantities, but this 
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Fig. 1. — Color plot of the density, logp in (gem -3 ) for the model with L = 0.6L e dd shown 
at different times given in years. Axes: horizontal: z: distance from equatorial plane in 
parsecs; R: distance from the BH in parsecs; 
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does not exclude the possibility of such coexistence on length scales smaller than 
we can resolve. 


5.1.2. Model with L = 0.5L e dd-' the importance of the hot flow pressure 

Figure 3 which shows a color plot of the density at different times shares a lot of similarity 
with Figure 1. One can distinguish approximately three regions in the flow: 1) a disk-like, 
geometrically thin and dense region, extending to R ~ 1.5 pc; 2) the IR-driven flow with 
lower density, separated by contact discontinuity from 3) a hot, photoionized flow. As L 
decreases, so does the energetics of the wind. However, the vertical thickness of the dense 
IR supported flow is determined by the delicate balance between a) the pv 2 pressure of the 
hot component from above, and b) by the pressure of IR supported wind from below. 

Decreasing T from 0.6 to 0.5 causes the pressure of the hot flow to be less than in the 
solution shown in Figure 3. The balance between the pressure of the photo-evaporated wind 
and IR drive dusty flow is is slightly shifted towards IR-dominated flow. 

Figure 4 shows the 2D distributions of the velocity components v z and vr. The slowly 
outflowing wind is clearly seen; its well defined boundary extends to 0.5 pc. The flow consists 
of a ’’core”, i.e. slow moving wind with radial velocities of a fewxlO < lOOkms -1 , and v z 
much smaller, of the order of the a fewxlkms -1 . Fast lower density component of the 
flow occupies most of the domain having velocities of the order of a few x 100 < 700kms _1 . 
The overall characteristics of the velocity field are similar to those of the model with L = 
0.6L e dd- Not surprisingly, pumping less energy initially in X-rays generate less powerful 
wind: Tkin/Aboi — lx 10 -5 — 7 x 10 -5 . 

The mass-loss rate (M) ~ 0.1M o yr _1 with the peak (M) ~ l.lM 0 yr _1 in the period 
of time between 8 x 10 3 yr and 10 x 10 3 yr, and the computation spans from 0 to 8 x 10 3 yr. 


5.1.3. Model with L = 0.3L e dd 

Figure 5 demonstrates the evolution of the density for this model at different times. The 
disk- wind system goes through several stages: the domination of the hot wind; puffing up 
the IR supported disk; squeezing of the later towards the equator by the hot wind; building 
a high density disk-like outflow, etc. As in previous examples, the system doesn’t come to a 
quasi-stationary state instead going through such episodes all over again. 

An interesting feature which is seen at t = 2.5 x 10 5 yr is in fact the residue from the 
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Fig. 2. — Color plot of the gas temperature, logT in K for the model with L = 0.6L e( jd, after 
4.2 • 10 5 yr. Axes: horizontal: z: distance from equatorial plane in parsecs; R: distance from 
the BH in parsecs; 

contact discontinuity which existed at earlier times. The average velocity (v) ~ a few x 
172kms _1 . This is about 60% of the escape velocity, [/ esc ; the numbers are given for t = 
20 t flyn . The maximum velocity of the dusty flow is 242kms _1 . 

The value of the mass-loss rate is remarkably similar to models with higher T : (M) ~ 
0.1 — O.2M 0 yr _1 reaching lM 0 yr _1 . 


5. 1. 4- Models with L < 0.1L e dd-' slow accretion with episodic outflow 

When the radiation input falls below an approximate threshold value no rigorous outflow 
is observed in our simulations. Instead we see that episodic outbursts of the hot evaporative 
flow are excited from the inner parts of the disk. Here we found this threshold luminosity to 
be L ~ 0.1L e dd- The results for the model for L = 0.1L e dd is shown in Figure 6. 

The density of this wind is too low to provide any considerable shielding from X-rays. 
Without much shielding it is difficult to form cold and dusty phase which can be accelerated 
by IR. Such an episodic outburst of the dense wind is followed by a gradual fall back of the 
dense component mostly in z direction towards equatorial disk. 
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Fig. 3. — Color plot of the density, log p in (gem -3 ) for the model with L = 0.5L e dd> shown 
at different times given in years. Axes: horizontal: z\ distance from equatorial plane in 
parsecs; R: distance from the BH in parsecs; 
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Fig. 4. — The surface plot of the 2 — and R— velocity components. Model is shown with 
L = 0.5L e dd, t = 6.5 x 10 5 yr; Horizontal: R: distance from the BH in parsecs; 2 : distance 
from the equatorial plane in parsecs; 
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The velocity of the hot wind, v < U esc , and it skirts the denser flow and falls back towards 
equator at larger radii, R ~ 2.5pc. The density piles up there and the shielding from X-rays 
rapidly increases. Correspondingly, the ionization parameter drops below £ m and the cold 
IR-driven wind again has a chance to develop. Notice that the vertical component of the 
radiation force is roughly g e ff, z ~ (dE /dr) cos 6*, where r is the spherical radius, and 9 is the 
inclination angle, measured from z axis. Since the vertical extent of the IR wind is small, 
geff^ is small, and the dynamic pressure of the IR wind cannot overwhelm the pressure of 
the hot component at higher z. The IR wind stalls and starts to slowly fall back (accrete) 
towards the equator. Density outside this dense IR-supported ’’pancake” drops and the gas 
gets overionized again. The cold disk itself shrinks from both inside and outside. 

We also calculated a model with L = 0.05Z/ e dd> and found that the general behavior of 
this model is similar to previous with L = 0.1L e dd- The wind is episodic and intermittent 
with the ” IR-supported pancake” phase. 

It is difficult to access the dynamical properties of models with L < 0.1L e dd> because 
our equatorial BC excludes the influx of matter into the accretion disk. Strictly speaking by 
assuming thin accretion disk as a source of matter for the wind together with outflowing BC 
conditions we assumed that the wind, even if very weak, is always present. We will further 
elaborate on this limitation in the following section. 


6. Column density and optical depth 

We calculate Thomson optical depth, tq = J' n e pdl , where l is measured along the line 
of sight at the angle, 9 from the vertical axis, and averaging is made over all calculated 
models for particular T. Processing different models, we are interested in an angle d p h where 
r p h = r(/9ph) = 1, i.e. when the wind becomes opaque. 

In the following we summarize results, listing models in the order from low to high 
luminosities. When T = 0.05 the optical depth rises sharply from 0.1 to 0.5 at (9) ~ 75°, 
(hereafter # e dge) and the photosphere is at (# p h) — 78°. 

Increasing the luminosity to T = 0.1 causes the dense part of the wind to become 
thicker in vertical extent, optical depth rises at the angle 9 ec j ge ~ 73°, and the photosphere 
is at (d p h) — 75°. At T = 0.3 optical depth rises sharply from 0.1 at (9) ~ 73°, (d p h) — 75°. 
and 9 e d g e — 73°. 

For the model with T = 0.5 the photosphere is at (9 P h) ~ 75° and the torus edge is at 
^edge — 73°. 
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Fig. 5. — Color plot of the density, log p in (g cm -3 ) for the model with L = 0.3L e dd shown 
at different times given in years. Axes: horizontal: z: distance from equatorial plane in 
parsecs; R: distance from the BH in parsecs; 
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Fig. 6. — Color plot of the density, logp in (gem -3 ) for the model with L = 0.1L e dd shown 
at different times given in years. Axes: horizontal: z: distance from equatorial plane in 
parsecs; R: distance from the BH in parsecs; 
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Figure 7 shows a color plot of a radial Thomson optical depth measured from a BH to 
a given point of the computational domain. The model shown has L = 0.6L e dd- This is the 
most luminous model we have calculated. It has the photosphere located at (9 P h) ~ 72° and 
the torus edge is at 6hd ge — 71°. 



Fig. 7. — Color plot of the radial Thomson optical depth, logr from a BH towards a given 
point. Axes: horizontal: z: distance from equatorial plane in parsecs; R: distance from the 
BH in parsecs; 


7. Mass-loss rate 

The opacity of dusty plasma in the infrared domain is 10-20 times that of 
electron scattering (Semenov et al. 2003). Such huge opacity makes possible for 
the AGN to have a strong wind, radiating at a fraction of Eddington luminosity 
(defined with respect to Thomson scattering (see Paper I for an estimate of an 
outflow onset) 

Mass-loss rate from the IR driven portion of the wind can be estimated 
combining the results from a theory of stellar winds and some of the approx- 
imations adopted in Paper I. From the stellar wind theory the mass-loss rate 
from a spherically symmetric wind can be found from the following approximate 
relation (Lamers & Cassinelli 1999): 
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M/(4vr) ~ L m /c (r IR - 1)/(4 tt u°°r IR )) r w , (14) 

where Ti R = Lx/L e dd,m? v°° is the wind terminal velocity, and r w is the wind optical 
depth in IR. 

For simplicity, as in Paper I in this section we assume that all incident 
UV and X-ray radiation is reprocessed within a narrow conversion layer which 
effective temperature, T e s can be found from the following approximate relation: 
ar^^edd = where T = L/L e jj, ry x = 0.5 is the fraction of X-ray radiation 

from the total radiation , and the fraction a ~ 0.5 of the incident flux is re- 
emitted outwardly. Adopting the above, one obtains Ti R = aTr/x^d/ K e which is 
approximately 1.25 for Mbh = 1O 7 M 0 , L = 0.5Lddd and K d = 10/c e . The value 
of v°° is estimated while neglecting P g in favor of the radiation pressure: v°° = 
y/2GM/r 0 (r IR — 1) ~ 207 km s -1 adopting the above set of parameters plus the wind 
launching radius R 0 = 0.5 pc, and assuming that the wind occupies a wedge-like 
(in a quarter of the domain) space of an opening angle of 75° we arrive to the 
estimate: M ~ 3 r w M 0 yr _1 . Monte-Carlo simulations show that the factor r w can 
reach 1-5 for an AGN radiating at a fraction of L e dd (Roth et al. 2012 ), meaning 
the potential for the wind tori to reach mass loss rate of 10M o yr _1 ; We do not 
see such mass-loss rates in our simulations, which may indicate that instead of 
being blown away as it would certainly do in a spherically-symmetric case, the 
gas is effectively ablated towards the equatorial plane by the pressure of the hot 
gas component. As the latter is itself the result of interaction with radiation, 
one can say that radiation effectively clears its way through the opening of the 
torus. 


8. Discussion 

The following average properties of our wind solutions can be summarized: 


1. The total mass- loss rate in the IR driven and warm absorber flows is between (M) ~ 
0.1M o yr _1 and (M) ~ 1.5M 0 yr _1 depending on the time of the evolution of the 
wind. This includes also the failed wind. The amount of gas which has enough energy 
to escape to infinity is in a good accord with the energetics derived from accretion. 
Notice that approximately 1.6 x 10 -1 — 1.3 x 10 -2 M 0 yr -1 are required to be accreted 
through a geometrically thin accretion disk to support our models. Also there is no 
problem of depletion of the torus since we have an infinite reservoir of matter in a razor 
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thin equatorial disk. 

2. The IR driven- wind has the velocity in the range of of a few x 10 — 200kms _1 , and the 
fast hot wind has velocity in the range of 400 — 800 km s -1 . 

3. The transmission properties are similar between different models. The characteristic 
angle at which the wind becomes opaque to Thomson scattering is about 72° — 75°, and 
is determined by the balance between the hot, photo- ionized wind and the cold wind 
supported by the IR pressure on dust. It is interesting that the geometry of the torus 
is close to one obtained in Dorodnitsyn et al. (2008) at late times when the pressure 
of the hot, evaporative component becomes important. 

4. Both hot and cold components of the wind are non-stationary. 

Models such as Konigl &; Kartje (1994); Elitzur & Shlosman (2006); Emmer- 
ing et al. (1992); Everett (2005); Keating et al. (2012); Bottorff et al. (2000) are 
developed on the premise of 1) the existence of a strong global magnetic field, 
and 2) uses the prescription of self-similarity. To a certain extent this limits the 
value of the direct comparison of our results with theirs. Indirect comparison 
would include: performing the detailed analysis of the spectral properties of the 
IR-driven ’’wind torus”; performing our simulations for the parameter range at- 
tributed to a well studied source (such as NGC5548, in case of Bottorff et al. 
(2000)). Those are the next things to do in order to assess the validity of the 
current model and both of these goals is a natural continuation of the line of 
present studies. We will elaborate on observational properties and the implications of 
the time-dependent behavior of our solutions in the following paper of this series. 


9. Conclusions 

Our studies of AGN infrared- driven winds have proceeded incrementally. In Paper I 
we considered a toy model which roughly approximated an infrared driven obscuring wind 
in AGN. An important prediction which resulted was that if the temperature of the dusty 
plasma in the rotating torus exceeds some characteristic value, the IR pressure on dust 
will inevitably overwhelm gravity, creating an outflow. However, this effective temperature, 
which approximately reads as T e g = — )^ 4 — 312 ( n ^ M? ) 1/ ' 4 — 987 (zmMx) 1 / 4 K, was derived 

analytically from a simplified model. 

In Paper I we combined calculations of stationary ID motion in the z-direction with 
2D calculations of the radiation field. The latter was calculated adopting 2D flux-limited 
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diffusion approximation. In Paper II we relaxed the assumption of ID stationary motion 
and adopted a 2.5D time-dependent picture. One important ’’toy-model” simplification still 
remained though: there was no X-ray radiation explicitly taken into account. Instead, we 
calculated the effective temperature of the gas at the boundary of the computational domain. 
In the present work we relaxed this assumption and included X-rays into the global picture 
explicitly. 

Despite these advances, our work still does not consider the following potentially im- 
portant physical processes: 

The assumption of a razor thin disk is idealized. An accretion disk at parsec scales is 
likely self gravitating, possibly affected by gravitothermal instability (i.e., Gammie 2001). 
This may provide a local source of energy and maintain a Toomre parameter 0 T ~ 1 Toomre 
(1964). In a marginally gravitationally unstable thin disk, the non-axisymmetric perturba- 
tions lead to angular momentum transport which can be described by an effective 
viscosity (i.e., Paczynski 1978; Shlosman et al. 1989, 1990; Rafikov 2009). Self- 
consistent models of geometrically thick disks in AGN at parsec scales call for the effects of 
self- gravity and radiation physics in 3D. This difficult endeavor can be a subject of future 
research. Another complication is that the problem of parsec scale winds may be intrinsi- 
cally connected with the problem of accretion at parsec scale. Relaxing the assumption of an 
equatorial thin disk and treating accretion adopting an effective a prescription is a possible 
next step in the development of the current model. 

The assumption of a razor thin disk as a source of matter underpins the current study. 
Certainly such a distribution of matter is in many respects most unfavorable for the formation 
of radiation driven flow, owing to the geometry of X-ray illumination. Thus, we expect that 
if the wind does form for certain T in the thin disk case, then we may expect that in less 
idealized circumstances it will also form. To relax the assumption of a razor thin disk one 
would need to abandon the assumption of equatorial symmetry and to include accretion disk 
into self-consistent consideration. 

Starting from the Paper I we see that much of the massive wind doesn’t leave the 
potential well of the BH. This conclusion is confirmed through Paper II, and in the present 
work as well. The precise fate of this gas is beyond the scope of these works, but the results 
suggest the following: if the accreted gas is rejected by the BH at least some of it will return 
for additional attempts. 

We conclude with a mix of results and ideas driven by our current studies: 

1. The distribution of radiation force inside the dusty component of the flow depends 
on the shape of the photosphere where X-rays are converted into IR and can only be 
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calculated in multi-dimensional simulations. This work includes both X-rays and IR 
into 2.5D time-dependent, hydrodynamic simulations . 

2. We do not observe a rotationally supported quasi-static torus in our simulations, al- 
though our simulated winds are likely to resemble some of the observed properties of 
such tori. 

3. The flow in the current simulations is more complex and time-dependent than that of 
our previous studies where we did not include X-rays explicitly. 

4. If geometrically thick dusty obscuration develops then a hot photo-ionized flow with 
velocities of 100 — 700kms _1 accompanies it. X-ray warm absorbers are the evapo- 
rative flow originating both from the disk and also evaporated from such cold dusty 
component. 

5. We speculate that large scale motions with v z = ±(a few) x lOOkms -1 originating 
from a spatially varying wind such we have calculated may be seen in maser emission 
in nearby bright type II AGNs. 
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Frequency-integrated moments E, F, which appear in the above set of equations (l)-(4) 
are obtained by calculating angular moments from the frequency-integrated specific intensity, 
/( r, 0, u, t ): 


10. Appendix A 



(15) 


( 16 ) 


This preprint was prepared with the A AS lAT^X macros v5.2. 
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The frequency-independent radiation pressure tensor P is found from 


P(r ,t) 



dQ nn/(r, 12, v, t ). 


(17) 


11. Appendix B 


Dust directly reprocess X-rays to IR and we approximately take this into account. The 
amount of energy absorbed by dust in a volume d 3 x: 


J_dEl 

d 3 x dt 


a n a 


faj^ddd^an^ c«., 


(18) 


where u x is the energy density of X-rays, a x = X^/ n <h and the number density of dust grains 
reads: rid = fdp/ m d, where m < \ is a mass of a grain, and / d is a dust-to-gas mass ratio. We 
make further simplifying assumption that the dust grain is being instantaneously heated to 
the effective temperature T f ] x . 

The contribution from the dust to the IR energy density (or, equivalently, to the tem- 
perature, T r is calculated from the following equation: 


—^ = n d ca d (T d;X - T r ) . (19) 

Assuming, that T d ,x is constant over the time interval [t, t + dt], the update for radiation is 
T djX (t + dt) = T d)X (l — e~ dt ! td ) -| — e~ dt ^ td T^ x {t). Notice, that E = aTf. The equation (11) 
is used to in operator-split fashion to update E prior to all other terms in equation (4) are 
updated. 


